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Abstract 

We study, via computer simulations, the fluctuations in the net electric 
charge, in a two dimensional one component plasma (OCP) with uniform 
background charge density —ep, in a region A inside a much larger overall 
neutral system. Setting e = 1 this is the same as the fluctuations in N\, the 
number of mobile particles of charge e. As expected the distribution of N\ 
has, for large A, a Gaussian form with a variance which grows only as /c|c?A|, 
where \dA\ is the length of the perimeter of A. The properties of this system 
depend only on the coupling parameter T = kT which is the same as the recip- 
rocal temperature in our units. Our simulations show that when the coupling 
parameter T increases, k(T) decreases to an asymptotic value k(oo) ~ k(2)/2 
which is equal (or very close) to that obtained for the corresponding variance 
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of particles on a rigid triangular lattice. Thus, for large T, the characteristic 



length £l = 2/c/p associated with charge fluctuations behaves very differently 



from that of the Debye length, £d ~ 1/y/T, which it approaches as T —* 0. 



The pair correlation function of the OCP is also studied. 
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I. INTRODUCTION 

A striking manifestation of the special long range nature of the Coulomb interactions 
and the resulting screening [[| is that the fluctuations of the net electrical charge Q\, 
contained in a subregion A of a spatially homogeneous and overall neutral equilibrium macro- 
scopic system, grow only as the surface area |<9A| and not as its volume |A|, the normal be- 
havior of fluctuations of extensive variables (at critical points the growth is even faster) ||- 
||. This behavior of the charge fluctuations can be readily understood by considering the 
(truncated) charge-charge correlation function S in an overall neutral translation invariant 
Coulomb system. This is given by, 



where p a (r) the microscopic density of particles of species a and e a their charge. The charge 
fluctuations of Qa is then expressed in terms of S as 




(1) 




(2) 



By introducing the characteristic function XA.( r ) of the domain A, 




0, r <£ A 



r e A 



(3) 



the integration in (FJ) can be extended to all space 
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where 



r)= [ dr' XA (r + r')[l- X A(r')} (5) 

JR d 



One observes now that the first term in (f|), which is proportional to the volume, vanishes 
when the integral of S vanishes, i.e. if there is "perfect screening" of the charges, 0- ]7). 
Under these conditions the only contribution to @ comes from the second term which yields, 
when the limit |A| — > M. d is taken in a self-similar way, 

<Ql> 



lim 



a d u d r d S{r)dr={-Y,P a el)i L (6) 
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|A|-+dc \dk\ 

where r = \r\, u d is the surface area of a unit sphere in ^-dimensions (u>i = 2), and it has 
been assumed that the system is also rotationally invariant. The geometrical constant a d 
is: «3 = 1/4, a 2 = 7r _1 , «i = 1/4, 0. The constant £l defines a characteristic length for 
the "charge separation" associated with the charge fluctuations [||]- [|J. 

We thus see that perfect screening and existence of the integral in (||) imply surface 
growth of the variance of charge fluctuations. The latter, but not the former, will generally 
be violated when a charged system with d-dimensional Coulomb interactions is confined to 
a lower dimensional space, in which case the fluctuations will grow as |<9A| log |A|, e.g. for an 
r _1 potential in d = 2 or logarithmic Coulomb potential in d = 1 @, [Kj. We also note that 
starting with a system which is not translation or rotation invariant, e.g. one with periodic 
structure, we can still apply (f|) and, when perfect screening holds, also (|) after averaging 
over translations and rotations. 

There are various arguments for expecting perfect screening, also referred to as the zeroth 
or charge "sum rule" , in equilibrium Coulomb systems - @ . They all involve assumptions 
about some minimal decay of correlations in such systems. This has been proven rigorously 
for classical systems at sufficiently high temperatures and low densities, i.e. in the Debye- 



Hiickel regime, where the decay is exponential, [[ll]], |I2]| . For quantum Coulomb systems the 
decay is only polynomial but still good enough. In fact one expects that perfect screening 
will always hold and so (Q\) /\ A| — ► as | A| — > oo. Of course in order to treat particles with 



charges of different signs via classical statistical mechanics it is necessary to modify their 
Coulomb potential at short range, e.g., by introducing hard cores to prevent collapse. This 
is so in dimensions d > 3 at all temperatures and in d = 2 at low temperatures. There is no 
such requirement for quantum systems as long as either the negative or positive charges (or 
both) obey Fermi statistics [|T^]. 

II. PARTICLE FLUCTUATIONS IN THE OCP 

A particularly interesting example of a system with reduced charge fluctuations is the 
one component plasma (OCP). In this system, used to model diverse physical situations, 
one kind of charges, say the negative ones, are treated as a uniform background with charge 
density —p, in which particles with positive unit point charges and average particle density p 
move about fL"4| , |L5|]. Since the OCP background is fixed, charge fluctuations correspond to 



fluctuations in particle number, N\. S(r) in (|J) now corresponds to the truncated particle- 
particle correlation function of a system with average particle density p i.e. S(r) = {pS(r) + 
p 2 [g(r) — 1]} in the conventional liquid theory notation. 

Surface area growth of the variance of particle fluctuations is a conceptually intriguing 
situation of interest beyond that of equilibrium Coulomb systems or even beyond statistical 
mechanics. Thus, it was proven in [T?J that for any system of point particles the variance of 



A^a, averaged over translations and rotations, grows at least as fast as \dA\ jT(|. Examples 
of systems having such a variance are particles arranged on a periodic lattice structure or 
those obtained via small distortions of such structures, [|16[], fll7| . It was actually proven re- 
cently that in one dimension such growth (which in d = 1 corresponds to bounded variance) 
implies, by itself, the existence of a periodic component in the extremal decomposition of 
any translation invariant measure, a fact already known for the d = 1 OCP, using directly 
methods of equilibrium statistical mechanics The existence of such a periodic compo- 
nent implies that there cannot be good decay of all correlations in a translation invariant 
state. In fact the only example (known to us) of fluctuations in a point particle system with 
good mixing (decay of correlation) properties in d > 1 is the OCP at high temperatures 



I^| . In d = 2 the OCP is exactly solvable when the coupling parameter T = (kT) 1 = 2 



T9fl . The truncated correlations between groups of particles separated by a distance D de- 
cay in this system in a super-exponential way, like exp[— cD 2 ], with c computed explicitly. 
(Interestingly the distribution of particles in the OCP at T = 2 is the same as that of the 
(suitably scaled) limiting distribution of eigenvalues of a random matrix M with entries 
Mij = Rij + ilij in which both and Iy are independent identically distributed Gaussian 
random variables 

In this note we study numerically the dependence of k on the coupling constant T, for a 
two-dimensional OCP. Using a unit of length proportional to p -1 / 2 the characteristic length 
£l = depends only on T. We expect on physical grounds that £l(F) will decrease with 
r so the question is: how small can the fluctuations become when T — > oo? On the one 
hand it is known from numerical studies that the 2d OCP undergoes some kind of ordering 



transition to a triangular lattice at T = T c ~ 140 [20], |2~1| . On the other hand the exact 
value of £z,(r) at T = 2 is only about twice the value it would have if the system was in 
a rigid triangular lattice and one averaged the fluctuations over translations and rotations 
[T7J . The question then is how will £l(P) behave as T increases towards and beyond r c ? 
For comparison the Debye-Huckel length which should approach as T — > ||, 

decreases as T -1 / 2 . If this, or something resembling it, was also the behavior of £l(P) we 
would have a system with fluctuations much below that of the rigid lattice, which would be 
surprising indeed. It is this question which motivated the investigations described here. 

The results of the simulations are given in Table I and Figure 1. They show unambigu- 
ously that the decrease in k(T) saturates for large T, approaching a value equal (certainly 
very close) to that of the triangular lattice. This is consistent with the intuition that there 
is a minimal value of the fluctuation per unit surface area and that this is achieved for a 
periodic arrangement of points. On the other hand it was recently shown that randomly 
distributing the position of each particle in Z d over a unit cell gives, for d ^ 350, smaller 
fluctuations than the rigid lattice P2"fl . The question of the minimal value of k in different 
dimensions is still open. 

We remark here that we have tried, so far unsuccessfully, to come up with a scheme 



for generating translation invariant, mixing measures of particle distributions in M. d which 
would have surface growth of the variance. These would be measures on points that do not 
start with an equilibrium distribution of the OCP having a finite number of particles in a 
box of volume V and then take the thermodynamic limit of V /* M d in an appropriate way 
to obtain an infinite particle system with average density p |TjJ . It is known that the infinite 
volume measure of the OCP is not Gibbsian because the probability of large deviations in the 
number of particles in a region A e M. d from its average value p|A| behaves like exp[— C|A| 7 ] 
with 7 > 1 p3| . One may wonder whether such behavior necessarily holds for all particle 
measures on M. d which have only surface growth of the variance. 
III. MODEL AND COMPUTATIONAL DETAILS 

In two dimensions the interaction between two particles of unit charge separated by a 
distance r is 

v(r) = -\n(r/L) (7) 

where L is an arbitrary unit length. A convenient unit of length, which will be used through- 
out the paper, is the radius of a disk containing one particle on the average, sometimes 
referred to as the "ion-disk radius", a ~ (np)~^. The reduced density is then p = n^ 1 and 
a thermodynamic state is uniquely defined by the coupling constant T. The difficulties as- 
sociated with computer simulations of this system due to the infinite range of the Coulomb 
interaction are well known. They are dealt with here, as in our previous work, by confining 



the particles to the surface of a sphere 



For N particles of unit charge moving on the surface of a sphere of radius R with uniform 



background density of opposite charge the total potential energy is taken to be |2T 

Vn = ln [zX^ 1 ~ Ui ' U ^\ ~ X l 1?\ (8) 

i<j 

where itj is a radial unit vector locating the position of particle i on the sphere surface. 
This corresponds to the distance between particles i and j being measured along the chord 
joining the particles. In the thermodynamic limit {N,R — > oo,p = N/AitR 2 constant) the 
energy differs from that of the planar system by a contribution of order 0(1/N) 



Most of our Monte Carlo simulations were performed with N = 1024 ions (i. e. R = 
y/N/2 = 16); some at low and high couplings used iV = 2048 (R = 22.62) to check the 
system size dependence. Charge fluctuations were calculated according to Eq. (^) which for 
the 2d OCP takes the form 

k = (< Nl > - < N A > 2 )/V A (9) 

where N A is the number of particles in the domain A drawn on the surface of the sphere and 
Va its perimeter. A convenient choice for A is a disk-like shape obtained by the intersection 
of the sphere with a cone of summit at the origin of the sphere and aperture 9. To check the 
independence of the results on surface shape additional computations were performed with 
a "rectangular" surface obtained by the intersection of the sphere with two planes parallel 
to and symmetric with respect to the equatorial plane of the simulation sphere and two 
parallel planes perpendicular to the equatorial plane. 
IV. RESULTS 
A. Charge Fluctuations 

Results for k in the range T = 0.01 — 140 covering the whole fluid domain are summarized 
in Table I and shown in Fig. 1. For each value of T, particle fluctuations were calculated for 
domains of different shape (disk- and rectangular-like) and different area. For instance, use 
of disk-like domains corresponding to 9 = 81.1°, 72°, 62.4° and 51.9°, and those symmetrical 
with respect to the center of the sphere, gave identical results (within statistical error) and 
were therefore averaged over. Results for rectangular domains (average over 8 domains) 
and disk-like domains were found to agree within statistical error (cf. Table I). When T is 
small k decreases rapidly but then saturates at a value k = 0.042 ± 0.002 near T = 80. For 
comparison, the corresponding value for k for particles on a rigid triangular lattice is (in our 
units) .0404 while for the square lattice it is .0411 fliTf . 

The excellent agreement of the simulation result with the exact value |18j k = 
(27T\/7r) _1 = 0.089793 at T = 2 (cf. Table I) shows that both system size and domain 
dimensions are sufficiently large for reliable results to be obtained. In fact system size de- 
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pendence is observed only for T^O.01. At these values of T both the 1024 and 2048 results 
differ from the Debye-Hiickel limiting value (27rv / 2T)~ 1 . This is presumably due to the fact 
that the correlation length (£i ~ £d) becomes comparable or exceeds the linear dimension 
of the domains. 

From Eq. it follows that k can be expressed in terms of the usual pair correlation 
function h(r) — g(r) — 1 

2 f°° 

k = — I dr r 2 h{r) (10) 

Jo 

where r is in units of a = (717?) -1 / 2 . 

An extremely good representation of k in the range < T < 2 can be obtained using in 
( PD an analytical approximation for the pair correlation function proposed in ref. . 

h{x) = -^{x^K^Xy/Ji) (11) 

leading to the simple expression 

1 1 r(/i + 3/2) 



K 



2n^^JI T(/i + l) 



(12) 



where fi = - — — and r(^) and K^{z) are the standard Gamma and Bessel functions, 



r 

respectively. It is exact at V = and T = 2 and reproduces the simulation data within 



statistical error. This is not so surprising since the zeroth, second and fourth moments of 
h(x) (Eq. (|il~l)) are exact (i.e. perfect screening, Stillinger-Lovett @ and compressibility 
sum rules are satisfied) and therefore an accurate value for the first moment Eq. (|i~l~D can 
also be expected. 

An analogous expression for T > 2 

1 1 r(z/) 

K ~ 277^7^ r>- 1/2) (13) 

based on Eq. (2.27) of ref. J24j]) {y = T/(T — 2) is less satisfactory. Although accurate 
near T = 2 the subsequent decrease is too slow giving (by analytic continuation) a value 
k = ^2 = 0.05066 in the zero temperature limit. 
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An accurate fit reproducing all the data within statistical error is given by 



K = 



0.11253954 1 + a^ 1 / 2 + a 2 r + a 3 r 3 / 2 

fva i + a 4 rV2 + r 



(14) 



(ai = 3.9896, a 2 = 1.1211, a 3 = 0.38138, = 4.0934). This form incorporates the exact 
Debye-Hiickel limit at low T and saturates at high values of T. 
B. Particle distribution 

As already alluded to in the Introduction it can be shown that under suitable clustering 
assumptions, which are fulfilled for the 2d OCP at small values of T (and at T — 2), the 
probability distribution of Qa/|<9A| _1 / 2 is Gaussian in the limit |A| — ► oo, i.e. the probability 
distribution P(N\) of particles in a domain A has Gaussian behavior ||, [|J. To this end a 
histogram of particle number N\ was recorded during the simulation and fitted to a Gaussian 
with variance a 2 . The distributions P(N\), calculated in a disk-like domain with 9 = 81.1°, 
and its Gaussian fits are shown in Fig. ^| for T = 0.5, 2, 10 and 100. We find, for example, 
after taking an average over different domains A, that a 2 /V = 0.165, 0.0901, 0.0551, 0.039 
for T = 0.5, 2,10,100, respectively, in good agreement with the direct calculation (cf. Table 



C. Screening lengths 

The Debye length £d measures the range of correlations between pairs of charges in 
the limit of high temperature ||, 0. For low values of F the asymptotic behavior of the 
pair correlation function exhibits exponential decay characterized by £d = l/v2T, 0], @- 
This length can be compared with the length £i typical of the spatial extension of charge 
fluctuations defined in (|6|). For the density fluctuations in our OCP has the form 



Table I shows that is close to £d for small values of T, the difference between £l and 
£d at T = 2 is ~ 10%. For T > 2 an estimate of the correlation length of the pair correlation 



I). 




(15) 
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function is more problematic. As shown in Fig. |^ the pair correlation functions have a 
damped oscillatory behavior at long range. These oscillations in h(r) can be represented 
within statistical error by the expression 

h(r) ~ Ae~ ar cos((3r - j)/y/r. (16) 

where A, a, 7 are fitting parameters. An example of such a fit is shown in Fig. § for T = 140. 
Eq. (|i~6|) can be obtained by expressing h(k), the Fourier transform of h(r) in terms of the 
direct correlation function c(k) according to h(k) = c(k)/(l — pc(k)) and assuming that the 
long range behavior of h(r) is driven by the poles in (1 — pc(/c)) _1 closest to the real axis 
("one" pole approximation studied, for instance, in |25| for the 3d OCP). The parameters a 
and 7 vary little with T. On the other hand a, which fixes the exponential damping of the 
correlations, decreases by a factor 3 between T = 20 and 140. Moreover, between T = 80 
and 140, a decreases by a factor 2 whereas k is nearly constant. It thus appears that for 
large V the correlations lengths £x, and £d are not related directly to the spatial decay of 
the envelope of h(r) in the OCP, c.f. ||. As mentioned above the exponential damping a is 
easily obtained from the "one pole approximation" but a more physical explanation of its 
origin has still to be found. 
V. CONCLUDING REMARKS 

Snapshots showing the arrangement of the ions on the sphere surface is given in Fig. |5| 
for T = 2 and 140. These show dramatically how the charges get more uniformly spaced as 
the temperature is lowered. On the other hand it would be hard to deduce from looking at 
the configuration at T = 2 (and even more so for smaller values of T) that the fluctuations 
were not normal: at T = 2, h(r) = — e~ 7rpr2 so (§) is certainly valid ||19|| . We therefore have 
to go beyond visual inspection and this can be done more easily for a multiple component 
system than for the OCP. 

As was noted in [Q surface area growth of the variance in Coulomb systems can be in- 
terpreted as corresponding to the tendency of charged particles to form "bound" neutral 
entities. This is most readily visualized in a two component system of charges ±e and densi- 
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ties p\ = P2 = p. Suppose now that these charges could be paired somehow to make neutral 
dipoles with bond length D. Then the charge fluctuation in a region A could be interpreted 
as resulting from the boundary OA, cutting across some of the "bonds" between the charges 
connecting the ions in these molecules. We would then get by standard arguments (assuming 
sufficient independence between molecules far apart) a central limit theorem for the fluctu- 
ations in Q\. The variance should then be as in @, i.e. of order CD e 2 p\dA\ with CD = 
The constant C measures the effect of correlation between the orientations of the dipoles 
cut by dA which is expected to be negative so C should be less than unity and decrease 
with T. The extension D might then be related to some screening length like £d at high 
temperatures and to a "hard core" or other minimal distance length at low temperatures 



and low densities 0, §, |6|. For the OCP D ~ a. 

It is of course not necessary in the above analysis to insist on pairs of charges, or dipoles, 
being the basic neutral entity. We could equally have neutral quadrupoles made up of two 
positive and two negative charges or some hierarchical structure as in the analysis of the 



Kosterlitz-Thouless phase in d = 2 || , [2E ] . What is necessary to make such an interpretation 
of charge fluctuations in A meaningful is that, in a typical configuration of the system, the 
neutral entities be spatially localized, i.e. not be greatly mixed up with other neutral entities. 
This is of course what happens in insulating materials, be they crystals, gases or liquids. 
Governed by quantum mechanics they consist of tightly bound neutral atoms or molecules. 
The picture in metals is similar in some ways to that of the OCP. Apparently enough of this 
picture remains true even in classical statistical mechanics of Coulomb systems, i.e. perfect 
screening, to give the correct behavior of the variance in Q\. 
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TABLES 

TABLE I. Particle number fluctuations k = (< > - < N A > 2 )/V in the 2d OCP as a 
function of L. V is the perimeter of the surface of A. One cycle consists of trial translations of 
the N ions, kdh = l/(27r\/2L), £d = l/y/2T and £l = 27tk. Length are in units of (7rp) -1 / 2 and 
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a the exact result is 0.089793 [HI], b results for "rectangular" surfaces. 
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TABLE II. Parameter a of the asymptotic decay of the pair correlation function 
h(x) ~ Ae~ ar cos(Ar — \pr (where r is in units of a the ion-disk radius). 
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0.278 


0.225 


0.194 
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FIGURES 

FIG. 1. Variation of k with coupling T 

FIG. 2. Particle number distribution P(N) for T = 0.05, 0.5, 2, 10 and 100 (from bottom to 
top). The filled circles are the Monte Carlo and the lines represent the Gaussian fits. 

FIG. 3. Pair distribution function g(r) for V = 10, 20, 40, 60, 80, 100, 120 and 140 (with 
increasing peak heights). 

FIG. 4. Fit of the long range oscillations of the pair correlation function h{r) = g(r) — 1 for 
r = 140 by means of the functional form Eq. 

FIG. 5. Snapshots of particle configurations on the sphere surface at T = 2, 140. Only the 
ions on the visible part of the surface are shown. (1024 particles were used for T = 2 and 2048 for 
r = 140.) 
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